Entanglement dynamics in a non-Markovian environment: an exactly solvable model 
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We study the non-Markovian effects on the dynamics of entanglement in an exactly-solvable model 
that involves two independent oscillators each coupled to its own stochastic noise source. First, using 
Lie algebraic and functional integral methods, we present an exact solution to the single-oscillator 
problem which includes an analytic expression for the density matrix and the complete statistics, 
i.e., the probability distribution functions for observables. We see non-monotonic evolution of the 
uncertainties in observables. Further, we extend this exact solution to the two-particle problem and 
calculate the entanglement in a subspace. We find the phenomena of 'sudden death' and 'rebirth' 
of entanglement. Interestingly, the time of death and rebirth is controlled by the amount of 'noisy' 
energy added into each single oscillator. If this energy increases above (decreases below) a threshold, 
we obtain sudden death (rebirth) of entanglement. 
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Noise in quantum systems can lead to abrupt and 
complete destruction (sudden death) of entanglement [l| . 
This is troubling since the death (sudden or otherwise) 
of entanglement represents one of the major obstacles 
towards building a quantum computer; see for example 
Q. In particular, when the bath is Markovian (memory- 
less as in [1]), the death of entanglement can be rather 
swift since the memory of the system's quantum state is 
wiped away by its totally uncorrelated interactions with 
the bath. 

Entanglement dynamics including sudden death and 
birth has been studied theoretically, e.g., in two qubits 
in several contexts [2, 0-01 an d m harmonic oscillators 
Q . The recent observation of these phenomena in pho- 
tonic systems Q and ensembles of atoms [§] has attracted 
great interest. In particular, it has been suspected that 
bath memory effects could not only provide an avenue to 
prolong entanglement but could also lead to its rebirth 
after it has experienced sudden death However, most 
noisy environments are hard to treat analytically by stan- 
dard techniques Q and one must use numerics or impose 
approximations to obtain a tractable result. 

In this Letter, we present an exactly solvable model 
involving two independent harmonic oscillators each in- 
teracting with its own classical non-Markovian stochas- 
tic reservoir. No back-reaction to the reservoirs is con- 
sidered. This system has the property that it can be 
solved analytically allowing us to study non-Markovian 
effects on the dynamics of entanglement including the 
prolonging of entanglement and its rebirth. Particularly, 
we study the dynamics of entanglement for the lowest 
two states of the oscillators which form a qubit-like sys- 
tem. Curiously, there is a one-to-one correspondence be- 
tween the amount of "noisy" energy added to each os- 
cillator and their entanglement: As the energy increases 
(decreases) across a threshold, we see sudden death (re- 
birth) of entanglement (see Fig. [1}. Furthermore, this 
initial-state dependent threshold is independent of the 
form of the noise correlations in time. 

The main advantage of our model is that it is exactly 
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FIG. 1. (Color online) Comparison between noisy energy 
of one of the oscillators and concurrence (entanglement) for 
baths with different memory. The initial state is (|01) + 
|10})/v / 2. When the energy exceeds (falls below) the thresh- 
old 0.455 a;, we see sudden death (rebirth) of entanglement. 
(E noise {t)) = ((E(t)}} ( - <(£(0)»s where ((E(t))) e is the ex- 
pectation value of the energy of a single oscillator averaged 
over noise. (a,b) use At = 0; (c,d) At — 4 and cj/A = 0.875; 
(e,f) At = 30 and cj/A = 0.25. 



solvable and non-Markovian 10]. We provide exact an- 
alytical expressions for a single oscillator's probability 
distribution function (PDF) of position, momentum, and 
energy where the memory of the bath is taken into ac- 
count explicitly. We find that the width of the PDFs of 
position, momentum, and energy are oscillatory in non- 
Markovian environments. Interestingly, this means that 
non-Markovian environments do not necessarily add un- 
certainty to the system monotonically. We also provide 
exact expressions for the correlation functions of observ- 
ables. In this situation we characterize the crossover from 
Markovian to non-Markovian behavior. 

Entanglement in harmonic oscillators can be quanti- 
fied in several ways and can be produced on demand 
with trapped ion systems 11| . Here, we focus on the low- 
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est two states of each oscillator which form a two-qubit 
Hilbert subspace. For a two-qubit system, entanglement 
is unambiguously quantified in terms of the concurrence 
C(p2(t)), where p 2 is the density matrix of two two-level 
system, we have 



C(t) = max{0, VAi 



A 4 }, (1) 



where Xi are the eigenvalues (in decreasing order) of the 
matrix p2(t)p~2(t) where p2_= { ~y® ~y)p2{ ~y® ~y)- Phys- 
ically, it can be shown [12J that states are maximally 
entangled if C(t) = 1 and completely unentangled for 
Cif) = 0. When C(t) = there exists a realization 
of p 2 (t) such that p 2 (t) = J2kPk IV'fc) (V^-l where every 
\ipk) is separable; i.e., the system is a classical mixture 
of separable states. The concurrence can vanish or ap- 
pear suddenly at a finite time, counter to what one may 
naively expect from the exponential decay of coherences 
(with characteristic time T%) which are local quantum 
phenomena. 

Consider the Hamiltonian of two harmonic oscillators 
H = Hi + i?2 in the presence of two statistically inde- 
pendent sources of classical noise (m = 1, 2) 



u{ala m + \) + ^ m \t)al + h 



(2) 



(no sum over m) where al n (a m ) are the standard cre- 
ation (annihilation) operators with [a n ,aj„] = S nm and 
£(™)(^) = £ { ™\t) +iti m) (t) are external stochastic fields 
with correlation functions of the form 



(^\t)6 m \t')) t = 5 nm K l3 (t,t') = S nm S tJ k(t - t') 



177,(5 - 
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where the average over noise is defined as the functional 
integral, ((• • • )) c = fV^ (- ; ■ ) exp[-| / /£ T if- l £] 
(suppressing the matrix indices for clarity). Our conclu- 
sions do not depend on the explicit form for the corre- 
lation function k(t — t') but on the time scale of decay 
of correlations t and the amplitude of the noise A. For 
t = 0, the bath has no memory and this leads to well 
known Markovian behavior [j| . We are mostly concerned 
with the regime where t ^ 0. 

We begin our analysis considering non-Markovian ef- 
fects on a single oscillator, i.e., Hi from Eq. ([2]). We cal- 
culate the PDF of position, momentum, and energy with 
Gaussian noise [10j . For an operator A its PDF P^(A\ t) 
is the noise averaged quantity -PJ^l; t] = (S[A— (A(t))])^. 
We find the expectation values of position and momen- 
tum ((i(t)) and (ft(t)) respectively) to be normally dis- 
tributed about the solutions of the classical equations 
of motion: X c \(t) = xq cos tut + pq sin ut and P c i(t) = 
— Xo sinwt + pq cos ujt. Given an initial coherent state of 
the form |V>(0)) = \zq) where zq = (xq-\-%pq)/v2, explicit 



calculation gives 
P*[X;t] = 
P p [P-t] = 



1 



^7rEx(*) 
1 



/ [X - X cl (t)] 2 \ 
CX P 1 ^ — ? . ( 4 ) 



exp 



2S.(t) 2 

[P - Pd(t)} 2 
2S p (t) 2 



(5) 



where X, P are random variables and the width of the 
PDFs is 

K, P (t)= f fdsds 1 g^&s^is, s')gf p (t,s'), (6) 
Jo Jo 

using Einstein summation convention from here on. Wc 
defined gf(t,s) = — s'muj(t — s), g p (t, s) — —cosuj(t~ 
s), g 2 (t, s) = cosw(i — s),g p (*) s ) = ~ sinuj(t — s). Notice 
that these Gaussian PDFs for both the position and mo- 
mentum are centered about the classical path and spread 
out in a Gaussian manner from that. From this, we can 
see that ((x(t)))^ and ((p{t)))^ will satisfy the standard 
classical equations of motion for the harmonic oscillator. 
Additionally, from Eq. ([6]) we can study the effects of the 
memory of the bath. For delta-correlated noise (r = 0) 
we find Brownian growth as Ti Xt p(t) — y/At which is ex- 
pected for Markovian-type of behavior. If we include 
memory in the bath we obtain 



ds ds' cos lj(s — s)k(s — s), (7) 
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which defines e(t), Fig.[2j this function appears frequently 
in our analysis. A change of variables gives 



du \u\k(u) cos urn. (8) 



e(t) = t J duk(u)cosu>u 
The first term in Eq. 



hence E 2 jP (i)) as t 



gives the slope of e(t) (and 
oo. For Gaussian time correlations 
k(t) = (A/tV^tt) exp[— t 2 /2r] we see that the slope is 
A exp[— w 2 r 2 /2] as t — > oo and the second derivative, at 
any time, is (2A/t^/2~k) cosuit exp[— i 2 /2r 2 ], so for short 
times the behavior is e(t) — > (2A/T\/27r)i 2 . Notice, the 
slope at long times vanishes exponentially as r increases. 

The non-monotonic behavior in S 2 p (i) means that 
the amount of spreading in the expectation value of x 
and p can decrease at short times if the bath is non- 
Markovian (r 7^ 0) counter to what one would naively 
expect given a noisy environment. At large times it 
grows linearly in time. This spreading is distinct, and 
in addition to the natural quantum mechanical spread 
of the system observables which is not captured by the 
PDFs. The non-Markovian effects can actually decrease 
the randomness in the system for short periods of time. 
In the extreme case of non-decaying noise correlations: 
Kij = ASij, £ 2 iP (£) = 2A(1 — cosLot)/uj 2 which means 
that the PDFs P$[X, 2-nkjuj\ and P p [P, 2nk/uj} (k is an 
integer) collapse into delta functions. At these discrete 
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times, the expectation value of these observables will 
yield the classical value of position and momentum with 
probability one and purely quantum mechanical behavior 
is restored. Intuitively, the system remembers its initial 
pure state and tries to restore it. When the memory is 
not long, the restoration is not complete but still can give 
non-monotonic behavior. 

We can also evaluate the PDF for the number operator 
analytically given an initial state Fock state \n). Explicit 
calculation [10( gives an expression that depends only on 



P A [N;t}= @iN - n) e-("-^\ 
e(t) 



(9) 



where N is the (continuous since it is expectation value of 
n) number random variable and O is the step function. 
Eq. ([9]) implies that only states with energy above the 
initial state \n) are occupied due to the heating. The 
exponential PDF has mean n + e(t) and variance e(t) 2 
with the memory of the bath being taken into account 
via e(t); see Fig.[5J The non-Markovian effects will cause 
the PDF to narrow as well, and in the limit of infinite 
time correlations of the bath it will periodically return 
to 5(N-n). 

All single oscillator statistical properties can be ob- 
tained from our PDFs. Non-Markovian effects also arise 
in the correlation functions of position and momentum. 
For example, if we let X(t) — [x(t), —p{t)] T and assum- 
ing that the initial state is \n) we can write Xi(t) = 
Xoi(t) + J* ds Rij(t, s)£j(s) where X (t) is the evolution 



without any driving force: xo(t) = (e 



andpo(£) = (e 



)/«V / 2. With these definitions 



the noise-averaged correlation functions are 

((X l (t)X J (t / ))) i = (X M (t)X Q] {t'))+ [ ds f ds' 

Jo Jo 

x R ^k {t, S )K kl ( S7S , )R l] ( S , ,t , ), (10) 

where Ru(t,t') = R 2 2(t,t') = cosuj(t - t'), Ri 2 (t,t') = 
— i?2i(i) = sincj(t — t'). (■ ■ • ) is the quantum mechanical 
expectation value and (•■■)? i s the average over noise. 
The average of the energy, ((E(t))) s /uj = (l/2)((x 2 (t) + 

p\t))) 6 = (1/2) Ei{{Ut)Mt))h ^ 
«£(t)» c /w - (n + 1/2) = (E noise {t))/uj 

= \ I [ dsds' tvR(t,s)K(s,s')R(s',t). (11) 
2 Jo Jo 

where we suppress matrix indices for clarity. Eq. (|11[) 
defines (E noisc (t)). Using Eq. © and defining A A = 
A — (A) we obtain 

((Apitf))^ = ((A.T(i) 2 )) ? = i(n+ 1/2) + e(t), (12) 




{E nois e(t)} = toe(t), 



(13) 



FIG. 2. (Color online) The function e(t) as defined in Eq. Q 
with Gaussian time-correlations. This quantity is a stand-in 
for: the rise in energy due to the stochastic forcing {E no ise(t)) 
(Eq. (|13l) ~). the variances in x and p averaged over time 
(Eqs. (1121) '). and the width of the PDF of position, momentum 
£^ p (t) (Eq. 0), and energy (Eq. ©). From top to bottom: 
cur = 3, ujt ~ 4.44, and lot = 6; in each case we have used 
u = A. 



We see linear Brownian-type behavior at large times (for 
all bath correlation times r) with the variance of posi- 
tion, momentum and excess (noisy) energy of the system. 
This is explicit in Fig. [IJi and Fig. [5] We also see oscilla- 
tions at short times for correlated baths [l3[ • Specifically, 
Fig. QJ shows an example for characteristic parameters 
ujt « 3.5 corresponding to moderate non-Markovian and 
Fig. with ujt w 7.5 corresponding to a strongly non- 
Markovian regime. We see a continuous cross-over from 
Markovian behavior, Fig. [Ij,, to strongly non-Markovian 
behavior, Fig. [T^, where the slope at long times vanishes 
exponentially, Fig. [21 

Next, we calculate the single oscillator density matrix 
to capture both quantum and statistical properties of 
our system. We start by exactly computing the density 
matrix p\ (t) for one oscillator averaged over noise realiza- 
tions. We define the density matrix using the quantum 
'Liouvillian' C(t) as 



h{t) = e £(t) pi(0) = {U{t)pi(p)W{t)) t , (14) 



where pi(0) = \ip(0)) (V'(O)I is t ne initial density matrix 
and U{t) is the single oscillator evolution operator which 
is the solution of id t U(t) = Hi(t)U(t), U(Q) = 1 and 
can be exactly solved [14|. The density matrix can be 
expanded as pi(0) = ^] nm a m „ \m) (n\ and therefore we 
only need to calculate the evolution of the basis elements 
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Explicit calculation [lCj gives 



. ^-~iu)t{rn—n) y 



e{t) m + k 2 F x [-k,- 



b~l+m,k+n 
k,l=0 

m; 1 + n — m, e(£) -2 



(n - m)!(l + e(*)) 



n+fc+l 



life! 



(15) 



for n > m and 2 -F\ is the hypergeometric function. For 
m > n exchange m ■<-> n and fc -f-> / on the right hand 
side of Eq. (fl~5j) (except in the e -* wt ( m ~ n ) phase factor) 
flp| . Next, we use this expression to compute the density 
matrix p(t) for two uncoupled entangled oscillators with 
Hamiltonian H = Hi + H 2 , Eq. © , 



p(o), 



(16) 



with the initial density matrix p(0) = |^o) (V'ol where 

|Vo) = ^(|01) + |10)), (17) 

and the states \nm) represents the first oscilla- 
tor in state |n) and the second in state |m). 
The density matrix can be written as p(t) = 
J2nm n'm' ( nm \p{t) \n'm') \nm) (n'm'\. We are only inter- 
ested in how the qubit-like entanglement in the subspace 
{ 1 00) , |01) , 1 10) , 1 1 1) } evolves in time. This defines anew 
4x4 density matrix p 2 given by p2(t) = Ilp(t)Il where 
II = Yin m=o \ nm ) (nm\ is the projection operator onto 
the subspace. We normalize this expression by the trace 
of il/5(i)n for convenience, but this does not effect our 
conclusions. Explicit calculation [l0| yields 
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(1 



e(t))3 

«t) 




V 





l+e(t) 

1/2 
l+e(t) 







1/2 
l+e(t) 

l+e(t) 









e(t)(l+e(t) 2 ) 

(l + £(t)P 



(18) 



Given this density matrix we compute the concurrence 
as given in Eq. ([1]). The results are presented in Fig. Q] 
along with plots of the energy added to a single oscil- 
lator. To see the connection between concurrence and 
energy, it can be shown directly that the energy given to 
a single oscillator by the stochastic forcing field is again 
Eq. (|T3|) (precisely: the energy is the average of the ener- 
gies for |0) and |1) time evolved separately). Eq. (fT%|) and 
Eq. (flU|) both depend on e(t), hence the energy controls 
the concurrence. Thus in Fig. Q] we show the behavior of 
the energy and concurrence on different time scales for 
the memory of the bath. We see that for a bath with 
no memory (white noise) the energy increasing linearly 
as a function of time and the concurrence vanishing at a 
critical time At c ps 0.455. For a bath with memory, non- 
Markovian oscillations of the energy in time lead to a re- 
birth of the entanglement as the energy crosses a specific 



initial-state dependent but noise-independent threshold 
(E no i se )/u> ~ 0.455. Intuitively, the system 'remembers' 
it quantum state, particularly its entanglement. This re- 
birth phenomenon is absent in baths with no memory 
(Fig. Ob). 

In terms of our initial density matrix, we may generate 
entanglement between higher energy states. Letting P — 
1 — II be the projection onto the rest of the Hilbert space, 
then schematically we have p = \ \p\ \ ■ I 'j>\ I • 1 \pl' ■ I'/il'. 
and only the first term Apil is separable when C(t) = 
(precisely: it can be written as the sum of density ma- 
trices of separable states) while the higher energy states 
may still exhibit entanglement between themselves and 
the lower energy states. Intuitively, the higher energy 
states act as a "cavity" to their respective "qubit" (as 
in [5|), so one may expect entanglement is being trans- 
ferred back and forth between them (as the classical noise 
slowly kills the total entanglement). 

In summary, we study analytically sudden death and 
rebirth of entanglement in a subspace of two entangled 
oscillators, and we showed that this is related to the 
memory of the bath. The longer the memory, the more 
the system 'remembers' its quantum state. Further, 
this effect is directly controlled by when the amount of 
"noisy" energy entering the system crosses a initial-state- 
dependent threshold. The same non-Markovian behavior 
is also responsible for a narrowing of the probability dis- 
tribution function of position, momentum, and energy of 
a single oscillator at short times. 
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SUPPLEMENTARY MATERIAL 



SYSTEM AND DEFINITIONS 



The Hamiltonian for a forced harmonic oscillator is (K = 1) 

H = w(ot + !) + [£(*)a f + T (*H 



(1) 



where £(t) = £i(t) + i&it) is a stochastic field, and a = -^(x + ip) where a\ a are the standard bosonic creation and 
annhilation operators. The noise correlator that fully characterizes our Gaussian noise is 



- K tJ (t,t') = 5 t] k(t - t 1 ). 
Averaging over the ^-fields can be done by path integration 

J V 2 i (• • • )e~5 & dt ' & dt"i 1 (t')K- j 1 {t',t")^(t") 



((• ••)>£ = 



fV 2 £e~^ % dv & *t"Si(t')Kr j 1 (*'<t")£i(.t") 



The evolution opertor for a single harmonic oscillator, id t U(t) — H(t)U(t) with U(0) = 1. is (14 



U(t) 



_ -iwt(ft+l/2)-i(*i(t)£+* 2 (t)£) t 7 (t) 



(2) 



(3) 



(4) 



where 3>i(i) = f Q dt' £j(t')Rji(t') (Einstein summation convention used throughout), h — a) a, and R is the 2x2 
rotation matrix 



R(t) = 



cos uit smut 
— sin cut cos cut 



(5) 



The phase 7(f) will be inconsequential for our calculations. It will further be useful define the adjoint of an operator 
ui jt Y=[X,t]. 



DENSITY MATRIX 



We compute the density matrix of a single oscillator by averaging over noise realizations of £i,2(£)- Using Eq. Q 
and e x Ye~ x — e ad *Y we obtain 



p(t;£(t)) = e -''" tft e —'( <I 'i( t ) i +*2(t)33) j g oe i(*l(t)i+<t'2(t)p) e «w*n 
_ g— iw* ad ft g— (t) ad £ +<t> 2 (t) ad p ) » 

where po is the initial density matrix. The noise-averaged density matrix is 



p(t) = (p(t;£(t))t 



-iwtadft / — i(<I>i (t) ad £ +<& 2 (t) adp) 



(6) 
(7) 



(8) 
(9) 
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Let x\ = x and X2 = p- Note further, [adf,adp] = and [ad#,l] = = [adp,l] allow us to treat ad^ and adp as 
c-numbers when integrating over £1.2: 



= g-iwtadfl J_ J V 2^ e ^P dt'Z % {t')R %] {t')^ ]e -±f t dt'p dt"U{t')K- ] \t'^ p o ( 1Q ) 

_ e -iut ad ft£; -i ad ii [J* dt' f* dt" ' Rij(-t')K jk (t' ,t")R kl (t")] ad i( ^ 

In Eq. (fTTj). we see the following operators showing up: {adri, ad? + ad?, ad? — ad|, 2 ada adp} which surprisingly, form 
a Lie algebra. This can be used to derive a full equation of motion for the density matrix. 
Schematically, Eq. (fTTj) can be written as 

Pit) = e £(t) Po, (12) 

where C(t) is the quantum 'Liouvillian' in some sense. Considering our particular form of noise in Eq. we can 
calculate 



Jo 



ad £i / dt' dt"R lJ (-t')K jk {t',t")R k i(t")&d £l =e(t)(adt + ad|), 



where 



e(t) = ds ds' cosa;(s — s')k(s — s'). 
Jo Jo 



(13) 



(14) 



aw m in 



and hence we need only consider how e c ^ acts on the state |m) (n|. First, we need a some 



facts about operators that act in this Hilbert space. It can be shown that any operator O can be expanded as 

dy dq 



O = 



2tt 



■ tiW 



yp-iqx] iqx-iyp 



Next, the operators e %qx lyp are eigenoperators of the operators ad^ and adp: 

ad £ e lqi - iyp = ye 1 ^-^, 
adp e lqi - lyp = qe lqii - ly P. 

And finally, we can calculate the matrix element 

(n\e^ q& \m) = ^(z*) m - n L^- n \\z\ 2 )e-^ 2 / 2 , z = ±(y + iq), 



(15) 



(16) 
(17) 



(18) 



where L n is an associated Laguerre polynomial. Also, it is easily shown that ad £ + adp commutes with ad^ (this 
can be seen from the Lie algebra). Combining these facts, we obtain from Eq. (ITBl that 



e cw \m) (n\ = e"^" 1 -") 



To evaluate this, we find the matrix element 

(k\{e cw \m) (n|}|Z) = e -^*( m -«) 

From Eq. ([18]), 



dy dq 
2vr 



' n\e iy v~ iqS: \m) e iq& ~ iyp e~ i e W( qS '+y 2 '> 



nW 
k\m\ 



dy dq 
2tt 



[n\e iyp - iq& \m) (k\e iq& - iy P\l) e-hW^+v 2 ' 



(19) 



(20) 



(21) 



Now, we let z = yjxe %e such that d 2 z — \dxdd 



(fc|{e £ « \m) (n\}\l) 



iW 



iutim-n) I dx I ^ x (fc-i+m-n)/2 el e(fc-;- m +r l ) i (m-«) (x)jL (fe-0 (a . )e -(l+ e (t)) a; ^2) 



(23) 



o 2^ 

S k -i, m - n e- lult{m - n) / dxx m ~ n L ( ,^ n \x)L ( i n ~ n) {x)e' {1+t{t))x . 
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using the identity 



to! n\ 



we can write 



(k\{e cw \m) (n\}\l) 



m\k\ 



5k-Lm- n e- tlJt{m - n) / dxx n - m L^- m Hx)L ( k n ' m \x)e- (1+e ^ x . 



(24) 



(25) 



The right hand side (RHS) of Eq. ([25]) is the same expression as the RHS of Eq. ([23| with n o to and k -o- I (except 
for the multiplicative e - lu)t (. m - n ) term). Thus, we can use Eq. (|23|) and assume m > n without loss of generality. At 
the end of our calculation, we simply switch indices to obtain m < n. 
A change of variables y = (1 + f (t))x yields 



(k\{e c(t) |m) (n\}\l) 



n\l\ 



Mm\ (1 + t{t)) m - n+1 
Togethcr with the property 



f dvv m ~ n L( m ~ n ) ( - I 



L 



(m— n) 
I 



IS 



l + e(t) 



L (:m-n) ( V \ _ 1 V- (t) n- 



i=0 



m \ T (m-n), , 



we obtain 



fZyy m - n L^ (- 



e(t) 



L 



(m—ri) 
I 



l + e(i) 



(1 


+ e(t)) n + l 




1 


(1 


+ e(t))"+' 




e(t) n+l 


(1 


+ e{t)) n + l 



m \ (I + m — n 
I- J 



EE^)" +i 

min(n,Z) , , 

e (t) n+l - 2l ( m ) ( + m ~ n \ ' ' 

i=o ^ n ^ ' ' 

ni \ (l-\-m—n 
I 



dyy m - n L { r n \y)L^- n \y)e-y 



(m — n)! 2 -P\ [— I, —n; 1 + to — n; e(t) 



where 2 F 1 is the hypergeometric function. Thus, we can return to Eq. (|26p to obtain 



(k\{e c ^\m) (n\}\l) = 



i\k\ & k -i„ eity^M-l-n-l + m-n-eity 



n\l\ {I + e{t)) m + l + l 
expanded in the number basis, we then have 

e £(t) |to) (n\ = e -^*( m -«) 



-iut(rn-n) 



(m — n)\ 



x < 



E 

1=0 

oo 

E 



\(l + m-n)\ e(t) n+l 2 F 1 [-l,-n;l 



i;e(ty 



(26) 



(27) 



(28) 

(29) 
(30) 

(31) 



n\l\ 



(m-n)!(l + e(t)) m + i + 1 



+ n - m)\ e(t) m+k 2 F 1 [-k, -to; l + n-m; e(i)" 



i!fc! 



(n-m)!(l + e(t)) 



n+fc+l 



+ to — n) (Z| for to > n, 



|fe) (fc + n — to| for n > m. 



(32) 



ENTANGLEMENT 



i independent oscillators with Hamiltonian 



Consider two 

= ujaW + ^ 1 >(t)a 1 +e i) (t)a\}+Loala 2 + (t)a 2 + ^ (t) 



! 4l 



(33) 
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i n) CO + ^2 n) (*) f or n = 1,2 and we impose (£<™W( m )^ 
|-0 O ) = ^7= (|0) |1) + |1) ® |0)) and the initial density matrix is 



where = ^ n) (t) + i^ n) (i) for n = 1,2 and we impose (C n) (t)d m) (t')) t = S nm 5ijk(t - t'). The intial state is 



p(0) = i (|0) (0| ® |1) (1| + |0) (1| ® |1) (0| + |1) (0| ® |0) (1| + |1) (1| ® |0) (0|) 



(34) 



Now assume that e^** 1 propagates the density matrix for the ith oscillator forward in time by t (as given by Eq. (|32p). 
so that p(t) is time-evolved by e £l ^ ® e £2 ^. We obtain 

/5(<) = i (e £l « |0> (0| ® e £2(t) |1) (1| + e £l ^ |0) (1| ® e £a(t) |1) (0| 

+e £lW |1) (0| ® e £2(t) |0) (1| + e £l W |1) (1| ® e £2(t) |0) (0|) . (35) 



From Eq. 



we have 



e AW| ) ( | = £ 



e(t)' 



(1 + £(*))'+ 



r 10 (i\ 



(=0 



(l + e(t))'+* 



e £ - ( * ) |0)(l| =e iwi ^VITl 



e(t)' 



2=0 



(1 +€(t))'+ 2 



; A(t) |x> (0| = e"^*^ vm 



6(t)« 



2=0 



(l + e(i))'+a 



10 + 



+ (Z| 



(36) 
(37) 
(38) 
(39) 



We are just interested in the entanglement between the bottom two states of each oscillator, hence if we call the 
projection operator tt = |0) (0| + |1) (1| and II = 7r ® tt, explicit calculation gives the 4x4 matrix 



flp(t)n 



(l + c(i)) 3 



(e{t) 


V 



o 

1/2 



l+e(t) l+e(t) 

1/2 H £ («) 2 

l+e(t) l+e(t) 

o g(«)[i+e(t) 2 



\ 



(40) 



(l+6(t)) = 



&(t) = n^n = 



The bath will mix all states, but Eq. (|40[) (properly normalized) is effecitively a two-"qubit" density matrix: 

iip(t)ii 

1 

) 2 + <t) + l) 

o \ 



tr np(i)n 



(4c(*) + l)(2e(i) 2 + c(*) + l) 

/e(i)(l + e(i)) 2 

(i +E ( t )2)(i + e(t)) |(l + e(t)) 



V 



|(l + e(t)) 




+ £ (t) 2 )(l + e(i)) 







(41) 
(42) 

(43) 



e(t)(l + eW 2 )/ 



Now, we can use p^ to calculate the concurrence by first computing p<i{£) — (a y ® a y )p2(t)(a y ® Cj,) and Q(t) = 
P2(t)p~2(t). The matrix Q has eigenvalues Ai(i) for i = 1,2,3,4 and where Ai(i) is the maximum eigenvalue. The 
concurrence is then 



C(t) =max{0, y/M(t)- v^-V^-V^} 



(44) 



The result is plotted in the main paper, Fig. 1. We see sudden death of entanglement whenever e(t) f=s 0.455090. 
Note that entanglement may persist with the higher states, but our calculation suggests entanglement between |0)®|1) 
and |1) ® |0) is lost. 
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PROBABILITY DISTRIBUTION FUNCTIONS 

In this section, we find the probability distribution functions (PDFs) of position, momentum, and energy. Consider 
an operator A(t), then its PDF is given by P^[^4;i] = (6 [A — (A(t))))^ where (A(t)) is the quantum expectation value 
of A. For example, given unspecified functions X(t) and gi(t), let 

(A(t))=X(t)+ f dt%(t') 9i (t'). (45) 



(This is the case for (x(t)} and (p(t)).) We then use the delta function identity 5(x) = J e lux (du / 2ir) to obtain 



P-\A;t] = [ — — [ T> 2 £(t)e~i & dtl ^ dt ^*(*i)^« 1 (*i.*2)£j-(t2)-/ t <Jt I «<(ti)ff < (ti)+iu(A-X(t)) < ^ 
J 2,7r Jv J 

The above is a quadratic path integral and can thus be solved exactly by standard techniques. K^ 1 represents the 
inverse integral kernel of K and we have used Kij(t,t') — Kji(t',t). Thus, 

P A [A;t} = j —e^^fo d 'i/o d*2fl < (ti)Jf«(ti,t2)fl J (t 2 )+t«(A-A'(t)) ^ 



1 



where the variance of this PDF is 



e -i(A-X(t)) 2 /£ 2 A (t) ( 4g ) 



Sx(*)= / dh f dt 2 g i (ti)K ij (t 1 ,t 2 )g j (t2). 
Jo Jo 



(49) 



The form of gi depends on A (see Eq. (l4"5"j) ). 
On the other hand, consider formally 



<B\I\) X(l)+ I dt 1 £ i (t 1 )g i (t 1 )+ [ dh f dt 2 ti(t 1 )Fi j {t 1 ,t 2 )Z j (t2). (50) 



o 



(this is the case for energy). In this case, 

P~\B;t] = [ ^ 1 / V 2 £{t)e~^ ti dtl /o dt 2^( t i)[^ 1 ( t i.*2)+2™^j(tl,t2)]? 3 (t2)-m/ t ciil«,(il)3i(il)+i«(S-X(t)) /g-g 

B J 2tx Vdetif 7 



2?r V det[l + 2iuifF] 



■ f* dt! f* dt 2 g t (ti)[K- 1 +2iuF]- j L (t 1 ,t 2 )gj(t2)+iu(B-X(t)) 



(52) 



The determinants can be viewed in the following way: To find "det D(t\, t 2 )" take the function D(ti,t 2 ) and time 
slice it iV — 1 times from to t, so one has an N x N matrix. Find the determinant of this matrix, then let N — > oo. 
1 + 2iuKF or K is also a 2 x 2 matrix, and in that case one just takes the deterimant of the matrix created by the 
direct product of those two spaces. This whole analysis hinges on the fact that Kij(tx,t 2 ) — Kji{t 2 ,tx) and further 
that Fijfah) = F ji (t 2 ,t 1 ). 

Position and momentum probablity distribution function 

If we let X% = x and X 2 = —p and we let the starting state be a coherent state \zq) with zq = -^(xq + ipo)i we 
have 

(Xi(t))=X iiCl {t)+ [ dt' Rijit-t')^'), (53) 
Jo 

where Xi lC \(t) are the solutions of the classical equations of motion for a harmonic oscillator, then from Eq. (j48|) we 
obtain 

Pa[ X-t] = ^J cxp(- [X ~ Xc 'y ]2 l , (54) 

x[ 1 V2^x x (t) I 2n x {ty /' y > 



1 expi [ P - P -'^ 2 
^V p {t) y \ 2E p (t) 



P P [P- 1] = - exp \ - l \^:Z n \ , (55) 
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and if we use Kij(t,t') — 5ijk(t — t'), we get S| (t) — e(t), where e(t) is given in Eq. (fT4|). 

Number (Energy) probability distribution function 

For what follows, we consider Kij(t,t') = Sijk(t — t') only. Given the initial state as a number state |n), we have 

(n(t)) =n + 

and using Eq. ([521) 



o Jo 



(h(t))=n + -l dh I dh&itjRijih-tzfefo), (56) 



iuKR] 

By functional integration methods, the determinant can be written as 



Now, let us consider the following quantity 

((0|{ e -«(*i(*) ad * +*2(«)ad j5 ) | j ( |}|0)) 4 = ( |{( e ^(*i(*)^ +»a(t)"d#}^ | ) ( |}| ) > ( 59 ) 

The left hand side of of Eq. (f59| can be found by standard techniques 

((0\{ e -^(t>iW^ +<Mt)a dj3 ) | j ( |}|0)) ? = ( e -^ 2 « Tfl «) c . (60) 

However, the RHS of Eq. (|5T)]) can be calculated just as in the Density Matrix section (see Eq. ©) with — > 2:$^) 
or equivalently e(i) — >• z 2 e(t). Reading off the |0) (0| component in Eq. we obtain 

(0|{(e^ ( * l( ' )ad *+* 2(t)ad ^) 4 |0)(0|}|0) =- l - T — (61) 

Letting z 1 = iu, we get 



1 1 



det [1 + mifi?] l + me(t)' 



(62) 



The PDF for the number operator is then 



P A [N;t]= I —- 7T = 4t/— TTT- (63) 

' 2tt1 + me(t) ie(t) J 2ir u - i/e(t) V ' 



Since e(t) > 0, this quantity is non-zero if N > n. Finally, we obtain 



P A [N;t] = e{N - n) e-^-^\ (64) 
e(t) 



as quoted in the main text. 



